Computational model for lipid binding regions in phospholipase (Ves a 1) from Vespa venom

The Thai banded tiger wasp (Vespa affinis) is a dangerous vespid species found in Southeast Asia, and its stings often result in fatalities due to the presence of lethal phospholipase A\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{1}$$\end{document}1, known as Vespapase or Ves a 1. Developing anti-venoms for Ves a 1 using chemical drugs, such as chemical drug guide, remains a challenging task. In this study, we screened 2056 drugs against the opening conformation of the venom using the ZINC 15 and e-Drug 3D databases. The binding free energy of the top five drug candidates complexed with Ves a 1 was calculated using 300-ns-MD trajectories. Our results revealed that voxilaprevir had a higher binding free energy at the catalytic sites than other drug candidates. Furthermore, the MD simulation results indicated that voxilaprevir formed stable conformations within the catalytic pocket. Consequently, voxilaprevir could act as a potent inhibitor, opening up avenues for the development of more effective anti-venom therapeutics for Ves a 1.

www.nature.com/scientificreports/ stings 13,14 . Anti-inflammatory drugs, antihistamines, and steroids are used to treat specific wasp sting symptoms, whereas anaphylaxis cases may require intrathecal epinephrine, intravenous epinephrine, and corticosteroids. However, caution should be exercised as hypertensive patients may be at risk of an acute hypertensive crisis, angina pectoris, or myocardial infarction 12 . Severe symptomatic treatments, such as respiratory support, cardiovascular support, haemodialysis, and haemoperfusion, may also be necessary. The hymenopteran insect, V. affinis, is known for its venom, which contains major protein allergens such as antigen 5, phospholipase A 1 B, hyaluronidase, and protease 5,15 . In 2013, the Ves a 1 protein was discovered in V. affinis venom 3 . The protein sequences of Ves a 1 and other related phospholipase A 1 protein families from different sources share a high degree of similarity with phospholipase A 1 (vPLA 2 ) protein from V. basalis and phospholipase A 1 VesT1.02 (VeT1.02) protein from V. tropica, with 93% and 92% similarity, respectively (Fig. S1B). All of these enzymes contain the typical catalytic triad residues Ser-X-Asn-X-His 4,16 (Fig. S1A,B).
The X-ray crystal structure of vPLA 2 shows that the α/β hydrolase fold at the N-terminus is necessary for catalytic activity. The tertiary structure dominants, β 5 and β 9, as well as the lid domain, are important in catalysis. The β 9 and lid domains play a role in substrate recognition and selectivity 4,17,18 . The catalytic triad residues of vPLA 2 are located between β 6 and α 8, β 7 and β 8, and α 9, respectively, as revealed by X-ray analysis. Ves a 1 also has an open conformation and a catalytic site similar to that of vPLA 2 from V. basalis 18 . Studies have shown that Ves a 1 has surface representation of catalytic triad (S170, D198, and H230) and auxiliary site (F241, Y242, N244, and Q249) 3,19 .
The activities of Ves a 1 are also involved in lipid membrane disruption 3,4 and haemolysis 20 . phospholipase A 1 s (PLA 1 s) family action binds the phospholipid bilayer and hydrolyses the sn-1 position of the phospholipid, eventually leading to cell lysis and the formation of lysophospholipids and free fatty acids 5,6,20,21 . Furthermore, recently proposed in silico modelling of phosphatidylcholines (PC) with venom PLA 1 s, with the short lid domain and β 9 initially interacting with the phospholipid polar heads. The hydrophobic residues in the β 5 immediately attract the substrate to the active site, which is formed by the S137-N165-H229 residues 22 . Interestingly, Ves a 1 was proposed as an auxiliary binding site to biological lipid membranes (Fig. 1), a conserved region for other hymenoptera venom PLA 1 s and mammalian lipases. The Ves a 1 auxiliary site could support Ves a 1 function since the protein was anchored on the phospholipid membrane surface with the help of a β 5 to promote its catalytic site improving the protein-lipid complex 3,17,22 . However, structural and molecular information about Ves a 1 venom remains unknown, and no specific anti-venom for Ves a 1 was found. As a result, our study attempted to identify critical atomistic information to be used for developing anti-venom for Ves a 1 in the state of chemical drugs remains a challenge for the study.

Methods
System preparation. The amino acid sequences of Ves a 1 was taken from Uniprot database (UniProtKB code: P0DMB4) and then submitted to AlphaFold 2 23 for structure prediction. Using the APBS web service 24 , the protonation states of all ionisable amino acids were assigned at pH 7.40. The ff14SB AMBER force field 25 and the GAFF2 (generalised AMBER force field) 26 . The TIP3P water 27  www.nature.com/scientificreports/ minimum distance of 10 around the protein surface and the solvation box edge, and the box dimensions were set to approximately 79.00 Å × 68.00 Å × 83.00 Å. The 29 sodium and 37 chloride ions were then added at random to neutralise at a physiological concentration of 0.15 M. The added hydrogen atoms and solvent molecules were subjected to 3000 steps of steepest descent (SD) energy minimisation with a total of 10,000 steps of minimisation to remove bad atomic contacts, while the rest of the molecules were held fixed. The protein and ligand were then minimised by 3000 SD over a total of 10,000 steps with constrained solvent molecules. Finally, the entire complex was minimised using the same method. The predicted free-ligand Ves a 1 was then simulated using molecular dynamics (MD) for 300 ns by AMBER16 28 16.00 Å × 20.00 Å × 13.00 Å for the auxiliary site was used, with the grid spacing at binding sites of each protein structure. The docking procedure was carried out using 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) (PubChem CID 5459377) as a control, which is a phospholipid commonly found in mammalian cell membrane 18 . The best predicted structure was used as a starting coordinate for MD simulation after competition.

Molecular dynamics (MD) simulation.
Prior to performing MD simulations, the partial atomic charges of the top five drug candidates were calculated using the restrained electrostatic potential (RESP) method at the HF/6-31G* level of theory with the Gaussian09 programmme 34 . Then, using the virtual screening, complexes of Ves a 1 with the control and the top five drug candidates were chosen as the initial conformation for MD simulation. The MD simulations were performed with the PMEMD.cuda module from the AMBER16 package 28,35 .
The system was equilibrated in a constant number (N), volume (V), and temperature (T) (NVT) ensemble after it was minimised to remove unusual inter-atomic contacts. The temperature in each simulated system was gradually increased from 0.01 to 310.00 K over 200 ps while a harmonic positional restraint of 10 kcal mol −1 Å −2 was applied to the C α atoms of protein. The particle mesh Ewald summation method was used to treat electrostatic interactions, with a cut-off distance of 9.00 Å for non-bonded interactions. All covalent bonds involving hydrogen atoms were constrained using the SHAKE algorithm 36 . Throughout the MD simulation, a simulation time step of 2 fs was used. The Langevin thermostat45 with a collision frequency of 1 ps-1 and the Berendsen barostat 37 with a pressure-relaxation time of 5 ps were used to control the temperature and pressure. To equilibrate the system, each complex was subjected to four steps of restrained MD simulations at 310.00 K with harmonic restraints of 5.00, 2.50, 1.00, and 0.10 kcal mol −1 Å −2 for a total of 800 ps, followed by another 200 ps with no restraint. Following that, each simulated system was run under the periodic boundary condition with the isothermal-isobaric of constant number (N), pressure (P), and temperature (T) (NPT) ensemble at 310.00 K and 1 atm pressure, similar to previous studies [38][39][40][41][42] until it reached 300 ns without restraint. The rootmean-square displacement (RMSD), solvent-accessible surface area (SASA), number of atomic contacts (#Atom contact), and number of ligand to Ves a 1 hydrogen bonds (#H-bond) were calculated from the number of Ves a 1 atoms within 3.50 Å of each drug. The distance between the hydrogen donor (HD) and acceptor (HA) of less than or equal to 3.50 Å and the HD-H · · · HA angles of more than or equal to 150.00 • were used to calculate the H-bond interactions. AMBER16's CPPTRAJ module and LigandScout software 43 were used to examine all trajectory analyses and pharmacophore interaction maps, respectively. The total binding free energy ( G SIE Bind ) was calculated using the solvated interaction energy (SIE) method, which was calculated as the following Eq. (1) 44,45 .
where E vdW is van der Waals, γ �SA is cavity, E Ele is electrostatic, and G RF is reaction field. E vdW and E Ele are denoted as intermolecular energies in the bound state. The coefficients in every calculation are α , γ , and C as 0.105, 0.013, and -2.89, respectively.
Free energy landscape (FEL) calculation. The free energy landscape (FEL) of Ves a 1 was calculated along with distance and the RMSD of T83-E93 and G251-C261 loops throughout the MD simulation. trial version where k b is the Boltzmann constant, T is the simulation temperature at 310.00 K, and g(x, y) is the normalised joint probability distribution. The minimum energy was set to zero. A bin size of 0.10 Å was used for the RMSD and distance. A total of 30,000 frames of 300 ns were used for the FEL calculations. www.nature.com/scientificreports/

Results and discussion
Structure and open/close conformation prediction of Ves a 1. The structure prediction was predicted by Alphafold 2 23 . The predicted structure from Alphafold 2 was then uploaded to MolProbity web services 46 to validate the predicted structures, which produced the Ramachandran plots, which were superimposed to vPLA 2 of V. basalis (PDB code: 4QNN), as shown in Figs. 1 and S2. The Ramachandran plot contained 96% of the preferred region, 100% of the allowed region, and 0% of the disallowed region (Fig. S2), and the RMSD of predicted Ves a1 compared to X-ray structure of vPLA 2 was 0.93 Å. The time evolution of the secondary structural elements of T83-E93 and G251-C261 loops (Fig. 2C right) showed that T83-E93 and G251-C261 loops were stable during simulation (300 ns, three-replicate), and root mean square fluctuation (RMSF) of T83-E93 and G251-C261 loops ( Fig. 2A,C right) showed high average values of fluctuations about 0.89 ± 0.23 Å, 0.96 ± 0.23 Å, and 0.87 ± 0.20 Å. This RMSF result may support the substrate recognition and substrate selectivity functions mentioned previously for G251-C261 loops, and this loop may act as a gate at the catalytic site as transitions between open and close conformations.
To gain a better understanding of the energetic changes associated with the conformational change of Ves a 1, the free energy landscape (FEL) was calculated as a function of distance and RMSD 47,48 of the T83-E93 and G251-C261 loops. The reaction coordinates of T83-E93 and G251-C261 loops were chosen as the reaction coordinates to monitor the structural change in order to obtain a two-dimensional energy landscape map of FEL. Because the T83-E93 and G251-C261 loops of Ves a 1 are superimposed, they are located at the same β 9 loop and lid-domain of vPLA 2 (PDB code: 4QNN). The FEL analysis results in Fig. 2D showed that the average distance values between the T83-E93 and G251-C261 loops of the two conformations were approximately 19.00 Å and 22.00 Å, indicating the close conformation and open conformation, respectively (Fig. 2B). As a result, the time evolution of the secondary structural elements (Fig. 2C left), RMSF, and FEL strongly suggested that the T83-E93 and G251-C261 loops could act as a gate at the catalytic site. The open conformation of Ves a 1 was then chosen as the initial structure for the following steps.

Results and discussion
Structure prediction and conformational analysis of Ves a 1. The prediction of the structure of Ves a 1 was performed using Alphafold 2 23 . The predicted structure was then validated using MolProbity web services 46 . The resulting Ramachandran plot revealed that 95.8% of the predicted structure was within the preferred region, 100% was in the allowed region, and 0% was in the disallowed region (see Fig. S2). The predicted Ves a 1 structure was superimposed with vPLA 2 of V. basalis (PDB code: 4QNN), resulting in an RMSD value of 0.93 Å(see Figs. 1, S2).
The stability of the T83-E93 and G251-C261 loops during simulation (300 ns, three-replicate) was confirmed by the time evolution of the secondary structural elements (see Fig. 2C right). The root mean square fluctuation (RMSF) of the T83-E93 and G251-C261 loops also indicated high average values of fluctuations, about 0.89 ± 0.23 Å, 0.96 ± 0.23 Å, and 0.87 ± 0.20 Å(see Fig. 2A,C right). These results support the substrate recognition and substrate selectivity functions of the G251-C261 loops, and suggest that this loop may act as a gate at the catalytic site by transitioning between open and closed conformations.
To further understand the conformational changes of Ves a 1, a free energy landscape (FEL) analysis was performed as a function of distance and RMSD 47,48 of the T83-E93 and G251-C261 loops. The T83-E93 and G251-C261 loops were chosen as the reaction coordinates to monitor the structural changes, and a two-dimensional energy landscape map of the FEL was obtained. The FEL analysis showed that the average distance values between the T83-E93 and G251-C261 loops of the closed and open conformations were approximately 19.00 Å and 22.00 Å, respectively (see Fig. 2B,D). This suggests that the T83-E93 and G251-C261 loops could act as a gate at the catalytic site, and supports the notion that the open conformation of Ves a 1 is the initial structure for subsequent steps.
Potential inhibitors of Ves a 1. The Ves a 1 protein was docked with 2,056 FDA-approved drugs and a ligand control (1,2-Dimyristoyl-sn-glycero-3-phosphocholine or DMPC) at both the catalytic and auxiliary sites using AutoDock Vina Version 1.2.0 32 (Fig. 1A). As the ligand control (DMPC) moved spontaneously from the auxiliary site to the catalytic site during all three individual simulations (Fig. 1B)

Conformational dynamics of Ves a 1 with DMPC and drugs complex.
The docking affinity scores of the top five Ves a 1 inhibitor candidates, including the control (DMPC), were evaluated using the AutoDock Vina Version 1.2.0. The ligand control and the drug candidates doxycycline, atovaquone, ubrogepant, dutasteride, and voxilaprevir interacted with Ves a 1 at the catalytic site, which was close to the solvent-exposed area. To evaluate the system's stability, the CPPTRAJ module was used to measure several parameters such as rootmean-square displacement (RMSD), solvent-accessible surface area (SASA), number of atomic contacts (#Atom contact), and number of hydrogen bonds (#H-bond) during the simulation (Fig. 4). All complexes showed high deviation in RMSD values for the first 30 ns and then fluctuated between 1.00 and 1.50 until the end of the simulation. The last 10 ns represented the equilibrium phase in each individual simulation (Run #1, Run #2, and Run #3). This phase was crucial for measuring the analysing of ligand control and the drug candidates. As the catalytic site of Ves a 1 is near the solvent-exposed area, water molecules may play a role in protein-ligand inter- www.nature.com/scientificreports/ actions. SASA calculations were performed on the residues within 5.00 Å of each ligand control and drug candidates to characterise the effect of water accessibility at the catalytic site of Ves a 1. The averaged SASA values were measured during the equilibrium phase of each individual simulation (Run#1, Run#2, and Run#3), the order of   www.nature.com/scientificreports/ SASA that referred to the region of catalytic site surface that is exposed sufficiently to interact with solvent molecules was DMPC (905. Binding affinity of of Ves a 1 with DMPC and drugs complex. The susceptibility of ligands (doxycycline, atovaquone, ubrogepant, dutasteride, and voxilaprevir) to Ves a 1 was estimated using the solvated interaction energy (SIE) method on 100 snapshots extracted from the last 10 ns of simulation. Table 1 shows the ranking of the binding free energy G SIE Bind (Run#1, Run#2, Run#3) from best to worst for voxilaprevir ( − 10.16 ± 0.43, − 9.81 ± 0.74, − 10.31 ± 0.40 kcal/mol) > ubrogepant ( − 7.48 ± 0.35, − 7.52 ± 0.32, − 7.32 ± 0.35 kcal/mol) > doxycycline ( − 6.97 ± 0.46, − 6.34 ± 0.41, − 7.03 ± 0.49 kcal/mol) > dutasteride ( − 6.29 ± 0.39, − 6.39 ± 0.62, − 7.07 ± 0.35 kcal/mol) > atovaquone ( − 5.89 ± 0.40, − 6.59 ± 0.32, − 5.88 ± 0.58 kcal/mol) and ligand control, DMCP ( − 9.49 ± 0.43, − 7.74 ± 0.64, − 9.34 ± 0.51 kcal/mol). The findings suggest that voxilaprevir may have a 1.10-fold (SIE) higher binding affinity to Ves a 1 than DMPC.
The MM/PBSA method was used to calculate G Residue Bind for key residues of voxilaprevir and DMPC, as shown in Fig. 5. Figure 5A,B depict the total energy contribution from each residue associated with ligand binding in both complexes, while Fig. 5C illustrates the drug orientation in the catalytic site of Ves a 1. Only residues with an energy stabilisation of lower than − 1.00 kcal/mol were considered. The results showed that six residues (F86, S88, L129, L171, P202, and L230) and eight residues (F86, S88, A92, L129, L230, I257, H263, and T264) were important for voxilaprevir and DMPC binding, respectively. Voxilaprevir and DMPC have the potential to interact with the catalytic pocket near the catalytic triad, S170, D198, and H263. The main energy contribution for stabilising voxilaprevir and DMPC could be electrostatic contribution in the range of approximately − 1.00 to 3.00 kcal/mol and vdW energy in the range of approximately less than − 2.00 kcal/mol, respectively.
The formation of hydrogen bonds between voxilaprevir and DMPC and their surrounding residues in the catalytic site of Ves a 1 may play a role in the inhibition of this targeted enzyme (Fig. 5C). www.nature.com/scientificreports/ www.nature.com/scientificreports/

Conclusion
After screening 2,056 drugs against Ves a 1, we selected the top five candidates for further investigation. We performed all-atom molecular dynamics simulations on these candidates for 300 ns. Then, we used SIE methods to estimate the binding free energy and understand the key binding interactions between voxilaprevir and Ves a 1. Our study identified the critical amino acid residues F86, S88, A92, L129, L230, H263, T264, and K268, which are involved in voxilaprevir's binding to the catalytic site of Ves a 1. The electrostatic contribution was found to drive the formation of the molecular complex between Ves a 1 and voxilaprevir. Additionally, hydrogen bond formations from pharmacophore models at S88, H263, T264, and K268 were identified as crucial for efficient binding. Therefore, voxilaprevir has the potential to act as a competitive inhibitor at the catalytic site of Ves a 1, opening up new possibilities for developing more effective anti-venom therapeutic agents for Ves a 1. To confirm the efficacy of voxilaprevir as an anti-venom therapeutic agent for Ves a 1, future experiments such as enzyme inhibitor assays and in vitro/in vivo studies may be necessary.

Data availability
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.